Observed2M_Model
Ben Weinstein
June 6, 2016
## [1] "Run Completed at 2016-07-29 14:53:50"#Load in data from Observed2m_Generate.Rmd
load("Observed.RData")
#load("ObservedModel.RData")
newModel<-T#source functions
source("Bayesian/BayesFunctions.R")View Raw Data
ggplot(indat,aes(x=Traitmatch,y=Yobs,col=as.factor(R))) + facet_wrap(~Hummingbird,ncol=4,scales="free") + geom_point() + geom_smooth(method = "glm",method.args=list(family="poisson")) + scale_color_manual("Resource Availability",values=c("black",'blue','red'))
ggplot(indat[indat$Yobs==1,],aes(x=Traitmatch,fill=as.factor(R))) + facet_wrap(~Hummingbird,ncol=4,scales="free") + geom_density(alpha=.7) + scale_fill_manual("Resource Availability",values=c("black",'blue','red'))
m<-glmer(data=indat[,],Yobs~Traitmatch*R+(1|Hummingbird),family="poisson")
summary(m)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: Yobs ~ Traitmatch * R + (1 | Hummingbird)
## Data: indat[, ]
##
## AIC BIC logLik deviance df.resid
## 6372.1 6414.8 -3179.0 6358.1 3314
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.0292 -0.6257 -0.3757 -0.2129 26.7497
##
## Random effects:
## Groups Name Variance Std.Dev.
## Hummingbird (Intercept) 1.536 1.24
## Number of obs: 3321, groups: Hummingbird, 19
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.610860 0.296018 -2.064 0.0391 *
## Traitmatch -0.032380 0.004328 -7.482 7.33e-14 ***
## RMedium -0.171362 0.094937 -1.805 0.0711 .
## RHigh -0.206707 0.088165 -2.345 0.0191 *
## Traitmatch:RMedium -0.003301 0.006738 -0.490 0.6242
## Traitmatch:RHigh -0.002527 0.006241 -0.405 0.6856
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Trtmtc RMedim RHigh Trt:RM
## Traitmatch -0.162
## RMedium -0.152 0.484
## RHigh -0.165 0.532 0.507
## Trtmtch:RMd 0.108 -0.646 -0.735 -0.347
## Trtmtch:RHg 0.119 -0.706 -0.347 -0.737 0.461
## convergence code: 0
## Model failed to converge with max|grad| = 0.00161504 (tol = 0.001, component 1)indat$pred<-predict(m,type="response")
ggplot(indat,aes(x=Traitmatch)) + geom_point(aes(y=Yobs)) + geom_line(aes(y=pred,col=as.factor(R)))
Hierarchical Occupancy Model
For hummingbird species i feeding on plant species j observed at time k and sampling event observed by transect (YTransect) or camera (YCamera)
Observation Model:
\[ YCamera_{i,j,k,d} \sim Binomial(N_{i,j,k},\omega_{Camera}) \] \[ \omega_{Camera} <- \phi_{Camera} * EffortCamera_k \]
Process Model:
\[ N_{i,j,k} \sim Poisson(\lambda_{i,j,k}) \] \[ log(\lambda_{i,j,k}) = \alpha_i + \beta_{1,i} * Traitmatch_{i,j} + \beta_{2,i} *Resources_k + \beta_{3,i} * Traitmatch_{i,j} * Resources_k \]
Priors
\[ \phi_{Camera} \sim U(0,1) \] \[\alpha_i \sim Normal(\alpha_\mu,\alpha_\tau)\] \[\beta_{1,i} \sim Normal(\mu_{\beta_1,\tau_{beta_1}})\] \[\beta_{2,i} \sim Normal(\mu_{\beta_2,\tau_{beta_2}})\] \[\beta_{3,i} \sim Normal(\mu_{\beta_3,\tau_{beta_3}})\]
Group Level Means \[ \mu_\alpha \sim Normal(0,0.0001)\] \[\mu_{\beta_1} \sim Normal(0,0.0001)\] \[\mu_{\beta_2} \sim Normal(0,0.0001)\] \[\mu_{\beta_3} \sim Normal(0,0.0001)\]
Group Level Variance \[\tau_{\alpha} \sim Uniform(0,1000)\] \[\tau_{\beta_1} \sim Uniform(0,1000)\] \[\tau_{\beta_2} \sim Uniform(0,1000)\] \[\tau_{\beta_3} \sim Uniform(0,1000)\]
#Source model
source("Bayesian/NmixturePoissonRagged2m.R")
#print model
writeLines(readLines("Bayesian/NmixturePoissonRagged2m.R"))##
## sink("Bayesian/NmixturePoissonRagged2m.jags")
##
## cat("
## model {
## #Compute true state for each pair of birds and plants
## for (i in 1:Birds){
## for (j in 1:Plants){
## for (k in 1:Times){
##
## #Process Model
## log(rho[i,j,k])<-alpha[i] + beta1[i] * Traitmatch[i,j] + beta2[i] * resources[i,j,k] + beta3[i] * resources[i,j,k] * Traitmatch[i,j]
##
## #True State
## S[i,j,k] ~ dpois(rho[i,j,k])
## }
## }
## }
##
## #Observation Model
## for (x in 1:Nobs){
##
## #Observation Process for cameras
## Yobs_camera[x] ~ dbin(dcam[Bird[x]],S[Bird[x],Plant[x],Time[x]])
##
## # #Assess Model Fit - Posterior predictive check
## #
## # #Fit discrepancy statistics
## #Camera
## eval_cam[x]<-dcam[Bird[x]]*S[Bird[x],Plant[x],Time[x]]
## E[x]<-pow((Yobs_camera[x]-eval_cam[x]),2)/(eval_cam[x]+0.5)
## #
## ynew_cam[x]~dbin(dcam[Bird[x]], S[Bird[x],Plant[x],Time[x]])
## E.new[x]<-pow((ynew_cam[x]-eval_cam[x]),2)/(eval_cam[x]+0.5)
## }
##
## #Priors
## #Observation model
## #Detect priors, logit transformed - Following lunn 2012 p85
##
## for(x in 1:Birds){
## #For Cameras
## logit(dcam[x])<-detect[x]
## detect[x] ~ dnorm(dprior_cam,tau_dcam)
## }
##
## #Detection group prior
## dprior_cam ~ dnorm(0,0.386)
##
## #Group effect detect camera
## tau_dcam ~ dt(0,1,1)I(0,)
## sigma_dcam<-pow(1/tau_dcam,.5)
##
## #Process Model
## #Species level priors
## for (i in 1:Birds){
##
## #Intercept
## #logit prior, then transform for plotting
## alpha[i] ~ dnorm(alpha_mu,alpha_tau)
##
## #Traits slope
## beta1[i] ~ dnorm(beta1_mu,beta1_tau)
##
## #Plant slope
## beta2[i] ~ dnorm(beta2_mu,beta2_tau)
##
## #Interaction slope
## beta3[i] ~ dnorm(beta3_mu,beta3_tau)
## }
##
## #Group process priors
##
## #Intercept
## alpha_mu ~ dnorm(0,0.386)
## alpha_tau ~ dt(0,1,1)I(0,)
## alpha_sigma<-pow(1/alpha_tau,0.5)
##
## #Trait
## beta1_mu~dnorm(0,0.386)
## beta1_tau ~ dt(0,1,1)I(0,)
## beta1_sigma<-pow(1/beta1_tau,0.5)
##
## #Resources
## beta2_mu~dnorm(0,0.386)
## beta2_tau ~ dt(0,1,1)I(0,)
## beta2_sigma<-pow(1/beta2_tau,0.5)
##
## #Interaction
## beta3_mu~dnorm(0,0.386)
## beta3_tau ~ dt(0,1,1)I(0,)
## beta3_sigma<-pow(1/beta3_tau,0.5)
##
## #derived posterior check
## fit<-sum(E[]) #Discrepancy for the observed data
## fitnew<-sum(E.new[])
##
## }
## ",fill=TRUE)
##
## sink()#Inits
InitStage <- function(){
#A blank Y matrix - all present
initY<-array(dim=c(Birds,Plants,Times),22)
initB<-rep(0.5,Birds)
list(S=initY,detect=initB)}
#Parameters to track
ParsStage <- c("alpha","beta1","beta2","beta3","alpha_mu","alpha_sigma","beta1_sigma","beta1_mu","beta2_mu","beta2_sigma","beta3_mu","beta3_sigma","dcam","dprior_cam","E","fit","fitnew")
#Jags Data
Dat<-list(
Yobs_camera = indat$Camera,
Birds=max(indat$jBird),
Bird=indat$jBird,
Plant=indat$jPlant,
Time=indat$jID,
Plants=max(indat$jPlant),
Times=max(indat$jID),
resources=resourceMatrix,
Nobs=nrow(indat),
Traitmatch=jTraitmatch)
#MCMC options
if(newModel){
system.time(
m2<-jags.parallel(data=Dat,parameters.to.save =ParsStage,inits=InitStage,model.file="Bayesian/NmixturePoissonRagged2m.jags",n.thin=4,n.iter=200000,n.burnin=200000-2000,n.chains=2,DIC=F)
)
}## user system elapsed
## 11.079 0.422 5267.675#recompile if needed
runs<-100000
recompile(m2)
if(!newModel){
system.time(m2<-update(m2,n.iter=runs,n.burnin=runs-2000,n.thin=2))
}#extract par to data.frame
pars_detect<-extract_par(m2,data=indat,Bird="jBird",Plant="jPlant")Assess Convergence
###Chains
ggplot(pars_detect[pars_detect$par %in% c("alpha","beta1","beta2","beta3"),],aes(x=Draw,y=estimate,col=as.factor(Chain))) + geom_line() + facet_grid(par~species,scale="free") + theme_bw() + labs(col="Chain") + ggtitle("Species Level")
ggplot(pars_detect[pars_detect$par %in% c("dcam"),],aes(x=Draw,y=estimate,col=as.factor(Chain))) + geom_line() + facet_grid(species~par,scale="free") + theme_bw() + labs(col="Chain") + ggtitle("Species Level")
ggplot(pars_detect[pars_detect$par %in% c("dprior_cam","beta1_mu","beta1_sigma","beta2_mu","beta2_sigma","beta3_mu","beta3_sigma","alpha_mu","alpha_sigma"),],aes(x=Draw,y=estimate,col=as.factor(Chain))) + geom_line() + theme_bw() + labs(col="Chain") + ggtitle("Group Level Parameters") + facet_wrap(~par,scales="free")
Posteriors
###Posterior Distributions
ggplot(pars_detect[pars_detect$par %in% c("alpha","beta1","beta2","beta3"),],aes(x=estimate)) + geom_histogram(position='identity') + facet_grid(species~par,scales="free") + theme_bw() + ggtitle("Species Parameters")
#Detection figure
pars_detect<-merge(pars_detect,jagsIndexBird,by.x="species",by.y="jBird",all=T)
ggplot(pars_detect[pars_detect$par %in% "dcam",],aes(x=par,y=estimate)) + geom_violin(fill='black') + theme_bw() + ggtitle("Detection Probability") + facet_wrap(~Hummingbird)
Detection Table
detecttable<-group_by(pars_detect,Hummingbird,par) %>% filter(par %in% c('dcam')) %>% summarize(mean=mean(estimate),lower=quantile(estimate,0.05),upper=quantile(estimate,0.95))
detecttable## Source: local data frame [19 x 5]
## Groups: Hummingbird [?]
##
## Hummingbird par mean lower upper
## (fctr) (fctr) (dbl) (dbl) (dbl)
## 1 Andean Emerald dcam 0.2756901 0.12795174 0.5023129
## 2 Booted Racket-tail dcam 0.2552547 0.16037808 0.3802274
## 3 Brown Inca dcam 0.2240020 0.12590289 0.3608807
## 4 Buff-tailed Coronet dcam 0.2284291 0.11936447 0.3819514
## 5 Collared Inca dcam 0.2559656 0.13416553 0.4183749
## 6 Crowned Woodnymph dcam 0.2782735 0.14273788 0.4410311
## 7 Fawn-breasted Brilliant dcam 0.2703968 0.10042306 0.4995664
## 8 Gorgeted Sunangel dcam 0.5285593 0.25151338 0.7797808
## 9 Green-crowned Brilliant dcam 0.2306932 0.07970214 0.4315912
## 10 Green-fronted Lancebill dcam 0.3600551 0.23995363 0.4987622
## 11 Hoary Puffleg dcam 0.3124265 0.15226076 0.5507799
## 12 Purple-bibbed Whitetip dcam 0.2735183 0.12942008 0.5228853
## 13 Rufous-tailed Hummingbird dcam 0.2671073 0.10285164 0.5073230
## 14 Speckled Hummingbird dcam 0.1924634 0.11030612 0.3612454
## 15 Stripe-throated Hermit dcam 0.2292337 0.14697289 0.3748659
## 16 Tawny-bellied Hermit dcam 0.2967347 0.18834135 0.4335916
## 17 Violet-tailed Sylph dcam 0.2612181 0.15820646 0.4405280
## 18 Wedge-billed Hummingbird dcam 0.1392272 0.05540372 0.2337274
## 19 White-whiskered Hermit dcam 0.2601769 0.17198555 0.3587414Discrepancy
The goodness of fit is a measured as chi-squared. The expected value for each day is the detection rate * the estimate intensity of interactions. The expected value is compared to the observed value of the actual data. In addition, a replicate dataset is generated from the posterior predicted intensity. Better fitting models will have lower discrepancy values and be Better fitting models are smaller values and closer to the 1:1 line. A perfect model would be 0 discrepancy. This is unrealsitic given the stochasticity in the sampling processes. Rather, its better to focus on relative discrepancy. In addition, a model with 0 discrepancy would likely be seriously overfit and have little to no predictive power.
fitstat<-pars_detect[pars_detect$par %in% c("fit","fitnew"),]
fitstat<-dcast(fitstat,Draw+Chain~par,value.var="estimate")
ymin<-round(min(fitstat$fit))
ymax<-round(max(fitstat$fit))
ab<-data.frame(x=0:ymax,y=0:ymax)
disc_obs<-ggplot(fitstat,aes(x=fit,y=fitnew)) + geom_point() + theme_bw() + labs(x="Discrepancy of observed data",y="Discrepancy of replicated data",col="Model") + ggtitle("Posterior predictive check") + geom_line(data=ab,aes(x=x,y=y)) + coord_fixed() + ylim(ymin=0,ymax=max(max(c(fitstat$fit,fitstat$fitnew)))) + xlim(xmin=0,xmax=max(max(c(fitstat$fit,fitstat$fitnew))))
disc_obs
#Bayesian p-value
sum(fitstat$fitnew>fitstat$fit)/nrow(fitstat)## [1] 0ggsave("Figures/ObservedDiscrepancy.jpeg",width = 5,height=10)Predicted Relationship
#Expand out pars
castdf<-dcast(pars_detect[pars_detect$par %in% c("beta1_mu","beta2_mu","beta3_mu","alpha_mu"),], Chain + Draw~par,value.var="estimate")Posterior prediction
#Trajectories from posterior
predy<-trajF(alpha=castdf$alpha_mu,beta1=castdf$beta1_mu,trait=indat$Traitmatch,resources=indat$scaledR,beta2=castdf$beta2_mu,beta3=castdf$beta3_mu)
ggplot(data=predy,aes(x=trait)) + geom_ribbon(aes(ymin=lower,ymax=upper),alpha=0.4,fill="red") + theme_bw() + ylab("Interactions") + xlab("Difference between Bill and Corolla Length") + geom_point(data=indat,aes(x=Traitmatch,y=Camera)) + geom_line(aes(y=mean)) + geom_point(data=indat,aes(x=Traitmatch,y=Transect)) 
At High and Low Resource Availability
#Trajectories from posterior
predH<-trajF(alpha=castdf$alpha_mu,beta1=castdf$beta1_mu,trait=indat[indat$R=="High","Traitmatch"],resources=indat[indat$R=="High","scaledR"],beta2=castdf$beta2_mu,beta3=castdf$beta3_mu)
predM<-trajF(alpha=castdf$alpha_mu,beta1=castdf$beta1_mu,trait=indat[indat$R=="Medium","Traitmatch"],resources=indat[indat$R=="Medium","scaledR"],beta2=castdf$beta2_mu,beta3=castdf$beta3_mu)
predL<-trajF(alpha=castdf$alpha_mu,beta1=castdf$beta1_mu,trait=indat[indat$R=="Low","Traitmatch"],resources=indat[indat$R=="Low","scaledR"],beta2=castdf$beta2_mu,beta3=castdf$beta3_mu)
predhl<-melt(list(High=predH,Medium=predM,Low=predL),id.vars=colnames(predH))
colnames(predhl)[5]<-"BFlowerL"
predhl$BFlowerL<-factor(as.character(predhl$BFlowerL))
levels(predhl$BFlowerL)<-c("Low","Medium","High")
ggplot(data=predhl,aes(x=trait)) + geom_ribbon(aes(ymin=lower,ymax=upper,fill=BFlowerL),alpha=0.5) + geom_line(aes(y=mean,col=BFlowerL),size=.8) + theme_bw() + ylab("Interactions") + xlab("Difference between Bill and Corolla Length") + geom_point(data=mindat,aes(x=Traitmatch,y=value))+ labs(fill="Resource Availability",col="Resource Availability") 
ggsave("Figures/AllRegression.jpeg",height=5,width=7)Species Predictions
castdf<-dcast(pars_detect[pars_detect$par %in% c("beta1","beta2","beta3","alpha"),], species +Chain +Draw ~par ,value.var="estimate")
#Turn to
castdf$species<-factor(castdf$species,levels=1:max(as.numeric(castdf$species)))
species.split<-split(castdf,list(castdf$species),drop = T)
species.traj<-list()
for(d in 1:length(species.split)){
x<-species.split[[d]]
index<-unique(x$species)
#get data for those species
billd<-indat[indat$jBird %in% index,]
#scale resources
species.traj[[d]]<-trajF(alpha=x$alpha,beta1=x$beta1,beta2=x$beta2,beta3=x$beta3,resources=billd$scaledR,trait=billd$Traitmatch)
}
names(species.traj)<-names(species.split)
species.traj<-melt(species.traj,id.var=colnames(species.traj[[1]]))
#split out names and model
species.traj[,c("Index")]<-colsplit(species.traj$L1,"\\.",c("Index"))
spe<-merge(species.traj,jagsIndexBird,by.x="Index",by.y="jBird")
#plot and compare to original data
ggplot(data=spe,aes(x=trait)) + geom_ribbon(aes(ymin=lower,ymax=upper),alpha=0.2,fill='red') + geom_line(aes(y=mean),size=.5) + theme_bw() + ylab("Daily Interaction Rate")+ xlab("Difference between Bill and Corolla Length") + facet_wrap(~Hummingbird,scales="free",ncol=4) + geom_point(data=mindat,aes(x=Traitmatch,y=value),size=2.5) 
Species Predictions: High and Low Availability
castdf<-dcast(pars_detect[pars_detect$par %in% c("beta1","beta2","beta3","alpha"),], species +Chain +Draw ~par ,value.var="estimate")
#Turn to
castdf$species<-factor(castdf$species,levels=1:max(as.numeric(castdf$species)))
species.split<-split(castdf,list(castdf$species),drop = T)
species.traj<-list()
for(d in 1:length(species.split)){
x<-species.split[[d]]
index<-unique(x$species)
#get data for those species
billd<-indat[indat$jBird %in% index,]
sh<-trajF(alpha=x$alpha,beta1=x$beta1,beta2=x$beta2,beta3=x$beta3,resources=billd[billd$R=="High","scaledR"],trait=billd[billd$R=="High","Traitmatch"])
sm<-trajF(alpha=x$alpha,beta1=x$beta1,beta2=x$beta2,beta3=x$beta3,resources=billd[billd$R=="Medium","scaledR"],trait=billd[billd$R=="Medium","Traitmatch"])
sl<-trajF(alpha=x$alpha,beta1=x$beta1,beta2=x$beta2,beta3=x$beta3,resources=billd[billd$R=="Low","scaledR"],trait=billd[billd$R=="Low","Traitmatch"])
sm<-melt(list(High=sh,Medium=sm,Low=sl),id.vars=colnames(sl))
colnames(sm)[5]<-"Resources"
species.traj[[d]]<-sm
}
names(species.traj)<-names(species.split)
species.traj<-melt(species.traj,id.var=colnames(species.traj[[1]]))
#split out names and model
species.traj[,c("Index")]<-colsplit(species.traj$L1,"\\.",c("Index"))
spe<-merge(species.traj,jagsIndexBird,by.x="Index",by.y="jBird")
#plot and compare to original data
indat$rplot<-as.factor(indat$scaledR)
levels(indat$rplot)<-c("Low","Medium","High")
ggplot(data=spe,aes(x=trait)) + geom_ribbon(aes(ymin=lower,ymax=upper,fill=Resources),alpha=0.2) + geom_line(aes(y=mean,col=Resources),size=.5) + theme_bw() + ylab("Daily Interaction Rate")+ xlab("Difference between Bill and Corolla Length") + facet_wrap(~Hummingbird,scales="free",ncol=4) + geom_point(data=indat,aes(x=Traitmatch,y=Yobs,col=rplot),size=4,alpha=.75)
ggsave("Figures/SpeciesRegression.jpeg",height=6,width=7)Species Level Interaction
castdf<-dcast(pars_detect[pars_detect$par %in% c("beta1","beta2","beta3","alpha"),], species +Chain + Draw~par,value.var="estimate")
#Turn to
castdf$species<-factor(castdf$species,levels=1:max(as.numeric(castdf$species)))
species.split<-split(castdf,list(castdf$species),drop = T)
species.traj<-list()
for(d in 1:length(species.split)){
dat<-species.split[[d]]
index<-unique(dat$species)
#get data for those species
billd<-mindat[mindat$jBird %in% index,]
#Calculate interaction effect
species.traj[[d]]<-intF(alpha=dat$alpha,beta1=dat$beta1,x=billd[billd$value > 0 & !is.na(billd$value),'Traitmatch'],resources=billd[billd$value > 0 & !is.na(billd$value),'scaledR'],beta2=dat$beta2,beta3=dat$beta3)
}
names(species.traj)<-names(species.split)
species.traj<-melt(species.traj,id.var=colnames(species.traj[[1]]))
#split out names and model
species.traj[,c("Index")]<-colsplit(species.traj$L1,"\\.",c("Index"))
spe<-merge(species.traj,jagsIndexBird,by.x="Index",by.y="jBird")
#match colnames
#plot and compare to original data
ggplot(data=spe,aes(x=x)) + geom_ribbon(aes(ymin=lower,ymax=upper,fill=Hummingbird),alpha=0.3) + geom_line(aes(y=mean,col=Hummingbird),size=1) + theme_bw() + xlab("Difference between Bill and Corolla Length") + ylab("Effect of Resources on Trait Difference") + facet_wrap(~Hummingbird,scales="free",ncol=3)
ggsave("Figures/SpeciesInteraction.jpeg",height=6,width=7)These plots can be tricky to interpret if one forgets that trait matching as a covariate is a distance. Therefore remember that a positive slope in the plot above indiciates, “As resources increase species use flowers less similiar to their bill lengths”.
Which species did we predict well?
smat<-pars_detect %>% filter(par=="E") %>% group_by(species) %>% summarize(fit=sum(estimate)) %>% arrange(desc(fit))
smat<-merge(smat,jagsIndexBird,by.x="species",by.y="jBird")
smat## species fit Hummingbird
## 1 1 4198.651 Andean Emerald
## 2 2 90919.540 Booted Racket-tail
## 3 3 209789.915 Brown Inca
## 4 4 38912.967 Buff-tailed Coronet
## 5 5 92914.052 Collared Inca
## 6 6 24295.666 Crowned Woodnymph
## 7 7 18216.200 Fawn-breasted Brilliant
## 8 8 31442.695 Gorgeted Sunangel
## 9 9 13621.877 Green-crowned Brilliant
## 10 10 41670.690 Green-fronted Lancebill
## 11 11 12204.590 Hoary Puffleg
## 12 12 12602.366 Purple-bibbed Whitetip
## 13 13 847.226 Rufous-tailed Hummingbird
## 14 14 95813.982 Speckled Hummingbird
## 15 15 94792.828 Stripe-throated Hermit
## 16 16 151151.034 Tawny-bellied Hermit
## 17 17 281231.597 Violet-tailed Sylph
## 18 18 33752.496 Wedge-billed Hummingbird
## 19 19 83972.069 White-whiskered HermitDiscrepancy matrix
dmat<-pars_detect %>% filter(par=="E") %>% group_by(species,plant) %>% summarize(fit=sum(estimate)) %>% arrange(desc(fit))
dmat<-merge(dmat,jagsIndexBird,by.x="species",by.y="jBird")
dmat<-merge(dmat,jagsIndexPlants,by.x="plant",by.y="jPlant")
head(dmat,10)## plant species fit Hummingbird
## 1 1 3 6192.83019 Brown Inca
## 2 1 2 124.05703 Booted Racket-tail
## 3 1 18 36.06558 Wedge-billed Hummingbird
## 4 1 9 50.77790 Green-crowned Brilliant
## 5 1 4 29.23103 Buff-tailed Coronet
## 6 1 10 259.95336 Green-fronted Lancebill
## 7 1 17 2121.24617 Violet-tailed Sylph
## 8 1 7 822.22259 Fawn-breasted Brilliant
## 9 1 16 1712.61751 Tawny-bellied Hermit
## 10 1 5 536.80132 Collared Inca
## Iplant_Double
## 1 Alloplectus tetragonoides
## 2 Alloplectus tetragonoides
## 3 Alloplectus tetragonoides
## 4 Alloplectus tetragonoides
## 5 Alloplectus tetragonoides
## 6 Alloplectus tetragonoides
## 7 Alloplectus tetragonoides
## 8 Alloplectus tetragonoides
## 9 Alloplectus tetragonoides
## 10 Alloplectus tetragonoides#get order
hord<-dmat %>% group_by(Hummingbird) %>% summarize(n=sum(fit)) %>% arrange(desc(n)) %>% .$Hummingbird
pord<-dmat %>% group_by(Iplant_Double) %>% summarize(n=sum(fit)) %>% arrange(n) %>% .$Iplant_Double
dmat$Hummingbird<-factor(dmat$Hummingbird,levels=hord)
dmat$Iplant_Double<-factor(dmat$Iplant_Double,levels=pord)
ggplot(dmat,aes(x=Hummingbird,Iplant_Double,fill=fit)) + geom_tile() + theme_bw() + scale_fill_continuous("Discrepancy",low="blue",high="red")
Resource Abundance Functions
Let’s take a closer look at distribution of interaction effect posteriors values for each species.
post<-pars_detect %>% filter(par %in% "beta2") %>% group_by(Hummingbird) %>% summarize(mean=mean(estimate),median=median(estimate),lower=quantile(probs=0.025,estimate),upper=quantile(probs=0.975,estimate)) %>% melt(id.vars='Hummingbird')
ggplot(pars_detect[pars_detect$par %in% "beta2",],aes(x=estimate)) + geom_histogram() + facet_wrap(~Hummingbird,scales='free',ncol=4) + geom_vline(data=post,aes(xintercept=value,col=variable))
Interaction density functions
Let’s take a closer look at distribution of interaction effect posteriors values for each species.
post<-pars_detect %>% filter(par %in% "beta3") %>% group_by(Hummingbird) %>% summarize(mean=mean(estimate),median=median(estimate),lower=quantile(probs=0.025,estimate),upper=quantile(probs=0.975,estimate)) %>% melt(id.vars='Hummingbird')
ggplot(pars_detect[pars_detect$par %in% "beta3",],aes(x=estimate)) + geom_histogram() + facet_wrap(~Hummingbird,scales='free',ncol=4) + geom_vline(data=post,aes(xintercept=value,col=variable))
Trait-matching and Bill Length
Do species with long bill lengths have positive traitmatching effects?
#species names
b<-pars_detect[pars_detect$par %in% "beta1",]
#traits
b<-merge(b,hum.morph,by.x="Hummingbird",by.y="English")
post<-b %>% filter(par %in% "beta1") %>% group_by(Hummingbird) %>% summarize(mean=mean(estimate),median=median(estimate),lower=quantile(probs=0.025,estimate),upper=quantile(probs=0.975,estimate),quantile_l=quantile(estimate)[[1]],quantile_u=quantile(estimate)[[2]]) %>% melt(id.vars='Hummingbird')
#get order of mean posterior
ord<-post %>% filter(variable=="mean") %>% arrange(value) %>% .$Hummingbird
b$Hummingbird<-factor(b$Hummingbird,levels=ord)
ggplot(b,aes(y=estimate,x=as.factor(Total_Culmen),group=Hummingbird)) + geom_violin(fill='grey50') + coord_flip() + ggtitle("Trait-matching and Bill Length") + theme_bw()
Interaction and Bill Length
Do species with long bill lengths have positive interaction effects?
#species names
b<-pars_detect[pars_detect$par %in% "beta3",]
#traits
b<-merge(b,hum.morph,by.x="Hummingbird",by.y="English")
post<-b %>% filter(par %in% "beta3") %>% group_by(Hummingbird) %>% summarize(mean=mean(estimate),median=median(estimate),lower=quantile(probs=0.025,estimate),upper=quantile(probs=0.975,estimate),quantile_l=quantile(estimate)[[1]],quantile_u=quantile(estimate)[[2]]) %>% melt(id.vars='Hummingbird')
#get order of mean posterior
ord<-post %>% filter(variable=="mean") %>% arrange(value) %>% .$Hummingbird
b$Hummingbird<-factor(b$Hummingbird,levels=ord)
ggplot(b,aes(y=estimate,x=Hummingbird,fill=Total_Culmen)) + geom_violin() + coord_flip() + scale_fill_continuous(low='blue',high='red') + ggtitle("Interaction Effect and Bill Length") + theme_bw()Estimated niche breadth
castdf<-dcast(pars_detect[pars_detect$par %in% c("beta1","beta2","beta3","alpha"),], species +Chain + Draw~par,value.var="estimate")
#Turn to
castdf$species<-factor(castdf$species,levels=1:max(as.numeric(castdf$species)))
species.split<-split(castdf,list(castdf$species),drop = T)
species.traj<-lapply(species.split,function(dat){
index<-unique(dat$species)
#get data for those species
billd<-indat[indat$jBird %in% index,]
d<-data.frame(alpha=dat$alpha,beta1=dat$beta1,beta2=dat$beta2,beta3=dat$beta3)
#fit regression for each input estimate
sampletraj<-list()
for (y in 1:nrow(d)){
v=exp(d$alpha[y] + d$beta1[y] * billd$Traitmatch + d$beta2[y] * billd$scaledR + d$beta3[y] * billd$Traitmatch*billd$scaledR)
sampletraj[[y]]<-data.frame(x=as.numeric(billd$Traitmatch),y=as.numeric(v),r=as.numeric(billd$scaledR),jBird=billd$jBird,jPlant=billd$jPlant,jTime=billd$jTime)
}
sample_all<-rbind_all(sampletraj)
})
species.traj<-rbind_all(species.traj)Mean Estimates for Corolla Sizes
species.mean<-species.traj %>% group_by(jBird,jPlant,r) %>% summarize(Traitmatch=unique(x),phi=mean(y))
tomerge<-indat %>% select(jBird,jPlant,Hummingbird,Iplant_Double) %>% distinct()
species.mean<-merge(species.mean,tomerge)
#get corolla sizes
species.mean<-merge(species.mean,fl.morph,by.x="Iplant_Double", by.y="Group.1")
#bill order
ord<-hum.morph %>% arrange(Total_Culmen) %>% .$English
species.mean$Hummingbird<-factor(species.mean$Hummingbird,levels=ord)
#add level to hum.morph to match naming convention
species.mean<-merge(species.mean,hum.morph[,c("English","Total_Culmen")],by.x="Hummingbird",by.y="English")
p<-ggplot(species.mean,aes(x=Corolla,y=phi,col=as.factor(r))) + geom_line(size=.9) + geom_vline(aes(xintercept=Total_Culmen),linetype='dashed') + facet_wrap(~Hummingbird,ncol=3,scales="free_y") + theme_bw() + ylab("Probability of Interaction") + scale_color_manual("Resource Availability",labels=c("Low","Medium","High"),values=c("Black","Blue","Red")) + xlab("Flower Corolla Length (mm)")
p
ggsave("Figures/ResponseCurves.jpeg",height=6.5,width=9)Niche Breadth
species.mean<-species.traj %>% group_by(jBird,jPlant,r) %>% summarize(Traitmatch=unique(x),phi=mean(y),phi_low=quantile(y,0.05),phi_high=quantile(y,0.95))
#merge names
species.mean<-merge(species.mean,jagsIndexBird)
species.mean<-merge(species.mean,jagsIndexPlants)
#get corolla sizes
species.mean<-merge(species.mean,fl.morph,by.x="Iplant_Double", by.y="Group.1")
#bill order
ord<-hum.morph %>% arrange(Total_Culmen) %>% .$English
species.mean$Hummingbird<-factor(species.mean$Hummingbird,levels=ord)
#add level to hum.morph to match naming convention
species.mean<-merge(species.mean,hum.morph[,c("English","Total_Culmen")],by.x="Hummingbird",by.y="English")
#label factor
species.mean$r<-as.factor(species.mean$r)
levels(species.mean$r)<-c("Low","Medium","High")
ggplot(species.mean) + geom_ribbon(alpha=0.7,aes(x=Corolla,ymin=phi_low,ymax=phi_high,fill=r)) + theme_bw() + facet_wrap(~Hummingbird,scales="free",ncol=4)+ ggtitle("Niche Breadth") + geom_vline(aes(xintercept=Total_Culmen),linetype='dashed') + geom_line(aes(x=Corolla,y=phi,fill=r)) + ylab("Probability of Interaction") + xlab("Corolla Length (mm)") + scale_fill_manual("Resource Availability",values=c("Grey70","Grey20","Black"))
ggsave("Figures/NicheBreadth.jpeg",height=8,width=9)Range plots
nsplit<-split(species.mean,species.mean$r)
makeR<-function(x){
#input matrix
aggm<-matrix(nrow=nrow(jagsIndexBird),ncol=nrow(jagsIndexPlants),data=0)
for (j in 1:nrow(x)){
aggm[x[j,"jBird"],x[j,"jPlant"]]<-rpois(1,lambda=x[j,"phi"])
}
aggm<-melt(aggm)
colnames(aggm)<-c("jBird","jPlant","P")
tomerge<-species.mean %>% select(jBird,jPlant,Corolla,Hummingbird,Traitmatch) %>% distinct()
aggm<-merge(aggm,tomerge)
return(aggm)
}
nstat<-lapply(nsplit,function(x){
netstat<-melt(lapply(1:50,function(k) makeR(x)),id.vars=c("jBird","jPlant","Hummingbird","Traitmatch","Corolla","P"))
colnames(netstat)[7]<-"Iteration"
return(netstat)
})
names(nstat)<-c("Low","Medium","High")
nstat<-melt(nstat,colnames(nstat[[1]]))Predicted density
ggplot(nstat[nstat$P==1,],aes(x=Corolla,fill=L1)) + geom_density(alpha=0.6) + facet_wrap(~Hummingbird,scales='free',nrow=5) + scale_fill_discrete("Resource Availability") + theme_bw()
rangestat<-nstat %>% filter(P==1) %>% group_by(Hummingbird,L1) %>% summarize(mean=mean(Corolla),lower=quantile(Corolla,0.05),upper=quantile(Corolla,0.95))
ggplot(rangestat,aes(x=Hummingbird,ymin=lower,ymax=upper,ymean=mean,col=L1)) + geom_linerange(alpha=0.5,position=position_dodge(width=.75),size=2) + geom_point(aes(y=mean),position=position_dodge(width=.75)) + labs(col="Resource Availability",y="Corolla Length (mm)") + theme_bw() + coord_flip() + theme(legend.position="bottom")
Predicted Mean
meanstat<-nstat %>% filter(P==1) %>% group_by(Hummingbird,L1,Iteration) %>% summarize(a=mean(Corolla))%>% summarize(mean=mean(a),lower=quantile(a,0.05),upper=quantile(a,0.95))
ggplot(meanstat,aes(x=Hummingbird,ymin=lower,ymax=upper,ymean=mean,col=L1)) + geom_linerange(alpha=0.5,position=position_dodge(width=.6),size=2) + geom_point(aes(y=mean),position=position_dodge(width=.6)) + labs(col="Resource Availability",y="Corolla Length (mm)") + theme_bw() + coord_flip() + theme(legend.position="bottom")
Predicted Min
meanstat<-nstat %>% filter(P==1) %>% group_by(Hummingbird,L1,Iteration) %>% summarize(a=min(Corolla))%>% summarize(mean=mean(a),lower=quantile(a,0.05),upper=quantile(a,0.95))
ggplot(meanstat,aes(x=Hummingbird,ymin=lower,ymax=upper,ymean=mean,col=L1)) + geom_linerange(alpha=0.5,position=position_dodge(width=.75),size=2) + geom_point(aes(y=mean),position=position_dodge(width=.75)) + labs(col="Resource Availability",y="Corolla Length (mm)") + theme_bw() + coord_flip() + theme(legend.position="bottom")
Predicted Max
meanstat<-nstat %>% filter(P==1) %>% group_by(Hummingbird,L1,Iteration) %>% summarize(a=min(Corolla))%>% summarize(mean=mean(a),lower=quantile(a,0.05),upper=quantile(a,0.95))
ggplot(meanstat,aes(x=Hummingbird,ymin=lower,ymax=upper,ymean=mean,col=L1)) + geom_linerange(alpha=0.5,position=position_dodge(width=.6),size=2) + geom_point(aes(y=mean),position=position_dodge(width=.6)) + labs(col="Resource Availability",y="Corolla Length (mm)") + theme_bw() + coord_flip() + theme(legend.position="bottom")
Generate network
#Split by resource
nsplit<-split(species.mean,species.mean$r)
makeN<-function(x){
#input matrix
aggm<-matrix(nrow=nrow(jagsIndexBird),ncol=nrow(jagsIndexPlants),data=0)
for (j in 1:nrow(x)){
aggm[x[j,"jBird"],x[j,"jPlant"]]<-rpois(1,lambda=x[j,"phi"])
}
#calculate network statistic
nstat<-networklevel(aggm,index=c("weighted connectance","weighted NODF","niche overlap"),level="lower")
}
nstat<-lapply(nsplit,function(x){
netstat<-melt(t(sapply(1:500,function(k) makeN(x))))
colnames(netstat)<-c("Iteration","Metric","value")
return(netstat)
})
names(nstat)<-c("Low","Medium","High")
nstat<-melt(nstat,colnames(nstat[[1]]))
ggplot(nstat,aes(x=value,fill=L1)) + geom_density(alpha=0.6) + facet_wrap(~Metric,scales='free',nrow=3) + scale_fill_discrete("Resource Availability") + theme_bw()
ggsave("Figures/NetworkStatistics.jpeg",height=5,width=6,dpi=600)Compared to raw visits
sindat<-split(indat,indat$scaledR)
ndat<-lapply(sindat,function(x){
web<-acast(x[x$Yobs==1,],Hummingbird~Iplant_Double,value.var = "Yobs", fun= function(x){ length(unique(x))})
nstat<-t(networklevel(web,index=c("weighted connectance","weighted NODF","niche overlap"),level="lower"))
})
names(ndat)<-c("Low","Medium","High")
ndat<-melt(ndat)
colnames(ndat)<-c("Var1","Metric","value","L1")
ggplot(nstat,aes(x=value,fill=L1)) + geom_density(alpha=0.6) + facet_wrap(~Metric,scales='free',nrow=2) + scale_fill_discrete("Resource Availability") + geom_vline(data=ndat,aes(xintercept=value,col=L1),linetype="dashed") + scale_color_discrete(guide='none')
save.image("ObservedModel.RData")